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We present a realistic model of the fragile glass former orthoterphenyl and the results of extensive 
molecular dynamics simulations in which we investigated its basic static and dynamic properties. 
In this model the internal molecular interactions between the three rigid phenyl rings are described 
by a set of force constants, including harmonic and anharmonic terms; the interactions among 
different molecules are described by Lennard- Jones site-site potentials. Self-diffusion properties 
are discussed in detail together with the temperature and momentum dependencies of the self- 
intermediate scattering function. The simulation data are compared with existing experimental 
results and with the main predictions of the Mode Coupling Theory. 

PACS numbers: 64.70.Pf, 71.15.D, 61.25.E, 61.20 



I. INTRODUCTION 

In the recent years a renewed interest on the glass transition phenomenon has motivated extensive experimental 
and theoretical works (see Q| and reference therein). On the theoretical side, new descriptions of the glass transition 
has been developed: they emphasized either the dynamic (as the Mode Coupling Theory (MCT) of Gotze or the 
coupled oscillators model of Ngai and Tsang Q|) or the thermodynamic (as the first principle computation based on 
a replica formulation of |^] or the inherent structure formalism computation of |||) aspects of the transition itself. 

A common feature of all these theories is they have been developed for "model systems" , often monoatomic models. 
The comparison of the theoretical results with the real experiment are, therefore, complicated by the trivial observation 
that in the real world the glass forming systems are made out of "molecular systems" . As a consequence in the current 
literature there is a large debate on the applicability of the theoretical predictions to the experimental outcome and 
the more stringent tests of the theories come from Molecular Dynamics (MD) works. 

As an example, it is highly debated in literature the origin, in molecular glassformer, of secondary relaxations that 
can be observed by several experimental techniques beside the well-known microscopic and structural dynamics (see, 
among others, [Q and reference therein). For instance, fast relaxations (i. e. in the 10 -12 s range) have been observed 
in several glasses and it is not yet clear if their origin is related to the molecular center of mass motion, as MCT 
would explain in terms of /3-process, or rather to rotational or intramolecular dynamics. It is clear the crucial role 
that MD simulations can play to solve this specific problem. If it is possible to build a realistic model able to take 
into account the internal degrees of freedom, as well as the translational dynamics, computer simulations allow one to 
access any observable quantity of the system, also those not directly measurable by present experimental techniques. 
Such possibility, together with the physical intuition, could allow the identification of the microscopic mechanisms 
underlying the different observed relaxation processes. 

In this paper we want to address the problem to set up a "realistic" potential for a glass model system capable to 
account for the internal molecular degrees of freedom. Among the glassforming molecular liquids characterized by an 
. ,_i \ extremely rich dynamical behavior, the organic fragile |IJ glass former ortho-terphenyl (OTP) (T m — 329 K, T c ~ 290 
K Q , Tq = 243 K) has received much attention from both experimental and numerical simulation points of view. The 
structure of OTP molecule, shown in Fig. [j], is known from neutron || and X-ray Jl(J diffraction studies; in condensed 
state the OTP molecules are bound together only by Van der Waals forces, which resemble the Lennard- Jones ones 
often used by most theories and computer simulations aiming to study the glass transition problem. 

Due to its structural complexity, if we would like to obtain reliable results to be compared with experimental data, 
we need to take into account not only the translations of the molecular center of masses and the rotations of the 
molecules as a whole, but also intramolecular motions like stretching along the molecular bonds, tilt of the bonds, 
rotation of the side rings with respect to the central one and so on. In other words, we need to describe the dynamics 
of the liquid at atomic level. On the contrary, in order to set up a computational scheme which is affordable in 
not-too-long time with the nowadays computer capability, we need the simplest model potential able to capture the 
relevant features of the dynamical behavior of the real system. 
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In the literature numerical studies of OTP have been proposed making use of different techniques ranging from 
harmonic lattice dynamics Jll| to molecular dynamics simulations on atomic level based on a general force field 
provided by the standard program Alchemy III p2|. Nevertheless, at our knowledge, only two studies based on 
molecular dynamics simulations of molecular models of OTP have been proposed so far: 

i) Lewis et al. jn^Jwj represent the molecule like a three-sites complex, each site playing the role of a whole phenyl 
ring, without internal dynamics and an intermolecular interaction of the Lennard- Jones (LJ) type. This model take 
in account only the dynamical behavior associated with the translations of the molecular center of masses and with 
the rotations of the molecules as a whole; 

ii) Kudchadkar and Wiest jl5| propose a more realistic model with the "true" structure of the molecule. The 
intermolecular interaction is of LJ type and, as internal degrees of freedom, only the rotational dynamics of the side 
rings with respect to the central one is taken into account. These internal degrees of freedom are effectively the most 
relevant, nevertheless the authors parameterize the potential in such a way that the side rings are, at equilibrium, in 
a configuration that corresponds to a saddle point in the molecular energy surface. 

The model potential we are going to introduce results to be much more efficient in mimicking the complexity of 
the dynamical behavior of the real system. The paper is organized as follows: in the next section we introduce the 
intramolecular model potential; in section ^ we explain how we calculated the force constants in order to reproduce a 
realistic isolated molecule vibrational spectrum. In section IV we present some computational details and in section |v| 
we remind some of the main well-established predictions of the ideal Mode Coupling Theory used in section [v| to test 
the center of mass dynamical behavior; In section VI we discuss our MD simulation results mainly with regard to the 
study of the diffusion and self-dynamic properties. The last section contains an overall discussion and the conclusions. 



II. THE MOLECULAR MODEL 



In our model the OTP molecule is constituted by three rigid hexagons (phenyl rings) of side L a — 0.139 nm 
connected as shown in Fig. |], i. e. two adjacent vertices of the parent (central) ring are bonded to one vertex of the 
two side rings by bonds whose length, at equilibrium, is Lj, = 0.15 nm . In our scheme, each vertex of the hexagons is 
thought to be occupied by a fictious atom of mass Mqh — 13 a.m.u. representing a carbon- hydrogen pair (C-H). The 
choice of such a fictious atom, with its renormalized mass, greatly simplify the computer simulation but presents some 
drawbacks: i) in the real molecule the two couples of carbon atoms connecting the rings are not bonded to hydrogen 
atoms, while in our model we consider all the 18 vertices having the same mass Mch so that the total molecular mass 
is overestimated (234 rather than 230 a.m.u) and ii) the moments of inertia of the model rings are smaller than the 
real ones the hydrogen atoms being too close to the ring center. Nevertheless, we expect only minor effects on the 
overall dynamics from the previous simplification. 

The three rings of a given molecule interact among themselves by an intramolecular potential, such potential being 
chosen i) to preserve the molecule from " dissociation" ; ii) to give the correct relative equilibrium positions for the 
three rings; Hi) to represent the real intramolecular vibrational spectrum as close as possible. The interaction among 
different molecules, actually among the rings pertaining to different molecules, is accounted for by a site-site pairwise 
additive potential energy of the (6-12) Lennard-Jones type, each site being one of the hexagons vertices. 

To sum up, the total interaction potential energy is written as the sum of an intermolecular and an intramolecular 
term: 

Vtot = Vinter + Vintra- (1) 

The first term can be written explicitly as: 

V mte r = - X! IZIZ v Lj(\rite - rjee'l) (2) 
ijij «' u> 

where f^i is the position of the ^-th atom (£=1. . . 6) in the £-th ring (£=1 ... 3, hereafter £=1 indicates the parent 
ring) belonging to the i-th molecule (i=l. . . N), and 



V LJ (R) = 4e 



r) ~ \r) 



(3) 



The optimal choice of the two intermolecular force parameters, e and a, will be discussed later. 

In principle, the intramolecular interaction potential can be expressed in terms of the degrees of freedom describing 
the center of mass positions (R^, £=2,3) and orientations (e.g. the set of Eulerian angles) of the side rings with 
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respect to the parent ring. However, for computational purposes, it is simpler to express the intramolecular potential 
in terms of orthonormal unit vectors attached to each ring or better in terms of quantities built from these vectors. 
With reference to Fig. ^ the sets of unit vectors for each ring j^, m^njj are defined as 

i 6 = (i,o,o) 

rh e = (0,1,0) 
n 6 = (0,0,1) 

i.e. 1% and rh^ are orthogonal unit vectors in the ring plane, while is the normal to that plane. 
The unit vectors that are parallel to the ring-ring bonds at equilibrium, are given by 




The positions of the four Carbon atoms that link the three rings, i.e Pi(2) an d Pi(3) m the parent ring and P 2 and 
P3 in the side rings, are given by, with respect to their ring centers R^, Pi (2) ~ Ri — -^a^i(2), Pi(3) — Ri — L a ui^, 
P 2 - R2 = L a u 2 and P 3 - R 3 = L a u 3 . 

Finally, it is useful to define four further interaction sites two pertaining to the parent ring, Si(2) = Ri + £ c Wi(2) and 
$1(3) = ^1 + L c Ui($) an( i two pertaining to the side rings (2 and 3 respectively), S 2 = R 2 — L c u 2 and S 3 = -R3 — L c u 3 . 
Here L c — L a + Lb/2, so that at the equilibrium position .§1(2) = S 2 and £1(3) = S3. 

The variable we have introduced will be used to simply express the different contributions to the intramolecular 
potential in the next subsections. 

A. Stretching along the ring-ring bonds and between the side rings 

The fluctuations of the distances among the centers of mass of the three rings are accounted for by introducing 
three "springs" . The parent-side ring stretching implies the elongation of a C-C bond, which is expected to have a 
high stiffness. The corresponding potential in harmonic approximation is written as: 

Vs = c 1 [\P m -P2\-L b } 2 + (4) 

ci [\p m - h\ - L b ] 2 

On the other hand, no direct chemical bond is present between the two side rings, and we expect a less stiff spring 
for the fluctuation of the distance between the side ring centers. We model this interaction by: 

V B =c 2 [\R 2 -R 3 \~{2L c )] 2 (5) 

The determination of the force constants c\ and C2, as well as the others we are going to introduce, will be discussed 
later. 



B. Tilt of the ring-ring bond 

In the OTP crystal structure the bond angles P 2 - Pi(2) - A(3) and P 3 - P 1{3) - P l{2) are 123.6° and 123.0° 
respectively, while the angles P 2 — A(2) — Pa an d P3 — A(3) — Pb are 118.4° and 117.4° (see Fig. |[) Further, in 
the isolated molecule, the ring-ring bonds are forced out of the plane of the parent ring so that the dihedral angle 
$ = P 2 — Pi(2) ~ Pi(3) ~ P3 1S 5.2°. This lack of planarity is due to the little asymmetry introduced by the difference 
between a Carbon bonded to an Hydrogen and a Carbon bonded to a Carbon of another ring. In our model, all these 
angles at the equilibrium are put equal to 120° and $ equal to 0. 
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We model the restoring forces for these angles by using the scalar product of the unit vectors u 2 and Ui(2) (as well 
as that of U3 and ui^). Being U2 ■ Ui(2) = "3 • ^1(3) = 1, at equilibrium, the quadratic term in the small oscillation 
approximation is given by: 

= c 3 (l - u 2 ■ ui(2)) + c 3 (l - u 3 ■ ui( 3) ) (6) 

However, this term is not enough to ensure the co-planarity of the vectors u 2 and U3 with the parent ring; to force 
the rings towards the co-planar equilibrium condition we make use of the "sites" S 2 and 5' 1 ( 2 ) (as well as S3 and §13) 
introducing between them a spring of vanishing equilibrium length: 

Vt 2 = c 4 \S 1{2} - S 2 \ 2 + (7) 

c 4|£l(3) — ^3 1 



C. Rotation of the side rings along the ring-ring bond 



Inside the intramolecular dynamics, we expect that a crucial role is played by the rotation of the side rings planes 
around the ring-ring bonds |l^,[l7j in interfering or modifying the intermolecular relaxation processes. The relevant 
variables describing this motion are the two angles {4> 2 , ^3} between the normals to the side ring and parent ring 
planes. In the crystalline structure, the two side rings have slightly different angles, 42.4° and 62.0° However we 
model the disordered condensed phases (liquid and glass) by using as equilibrium angle values those of the isolated 
molecule, i.e. 54.0°, as discussed below. We have to remark that the isolated molecule symmetry implies two iso- 
energetic configurations separated by a finite barrier as is qualitatively illustrated in Fig. ||. We performed an ab-initio 
calculation of the single molecule potential energy surface as a function of (f> 2 and ^3 with all the other internal degrees 
of freedom fixed to their equilibrium values. Such calculation consists in the minimization of the Hartree-Fock energy 
over the gaussian basis set 3-21G. For each atomic species the inner shell is made up of 3 gaussians while the valence 
shell is a linear combination of 2 gaussians orbital plus a 1 gaussian orbital; the minimization was carried out by the 
standard package Gaussian 94. In Fig. || a contour plot of the Hartree-Fock potential energy is shown as a function 
of the rotation angles. It is important to quote that this map has only a semi-quantitative meaning, the structure of 
the whole molecule being not re-optimized during the scan of 4> 2 and <p3 angles. However, a careful study, performed 
reoptimizing the whole molecular structure, has been done around the saddle point (f> 2 — <p3 — 90° and around the 
two equivalent minima that turn out to be at 4> 2 — 4>3 — 54° and 4> 2 — 4>3 = 126°. From this calculation the barrier 
height has been estimated to be V s /k,B = 580 K and, from this value, it is possible to envisage the nature of the 
rotational motion at the temperatures of interest: the two side rings can pivot in phase around the bonds crossing 
from one minimum to the other degenerate one. Moreover they can perform librational out-of-phase motions of 
(approximatively) harmonic type. 

In order to represent this potential surface we express the in-phase rotation of the two side rings with a high-order 
(6th) polynomium and the out-of-phase rotation in the harmonic approximation. For this purpose we use as primary 
variables the scalar products ri\ ■ h 2 and hi ■ ¥13 (their equilibrium value being (hi ■ h 2 ) eq = (hi ■ h 2 ) eq = uq = 0.59) 
and to disentangle the in-phase from the out-of-phase motion, we introduce the two variables: 

hi ■ h 2 + hi ■ h 3 

a = 2 (8) 

m ■ n 2 - m • n 3 

P = 2 

The final form of the " internal rotation" potential will be 

V R = V Rl (a)+V R M (9) 

with 

V Rl (a) = ha 2 + b 2 a 4 + b 3 a e (10) 

V R2 (/3) = c 6 p 2 (11) 

The parameters 61, b 2 and 63 describing the in-phase rotation potential are derived according to the following 
procedure. In the harmonic approximation, in proximity of the equilibrium position a Q for the scalar products, it 
must hold 
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V Rl (a)~c b {a-a ) 2 ; (12) 

and taking into account that the barrier height must be equal to V s , we have the following conditions: 

V Rl (a ) = -V s 

V Rl (a o )=0 

VH 1 (a )=2c 5 



implying 



& 2 = 77^-3W (13) 



b. 
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III. FORCE CONSTANTS 



We can finally write our internal model potential like 

V mtra = V S + V B + V Tl + V T2 + V Rl + V R2 ; (14) 

this potential is parametrized by the set of six free coefficients {ck} whose actual value can be tuned in order to obtain 
a realistic free molecule vibrational spectrum. Indeed, in the small oscillations approximation, we can determine the 
values of the coefficients {c^} by diagonalizing the dynamical matrix and fitting the resulting eigenfrequencies Up IAG 
to the lowest frequencies uj^j F obtained by an Hartree-Fock calculation of the vibrational frequencies in the electronic 
ground state of the isolated molecule. 

We have 18 eigenvalues for the dynamics of three rings (A = 1, ... 18) but actually only iV e = 12 eigenvalues uip IAG 
are non vanishing, the other being associated with the translations rotations of the molecule as a whole. 

More explicitly, the set {c^} is obtained by minimizing, by the standard Levenberg-Marquardt algorithm fTs|] , the 
error function 

N e 

Y^^DIAG-^HF? (15) 
A=l 

where, for a given set of {c&}, the quantities u>p IAG are the solution of the eigenvalues problem Ji"9| 

\V-lu 2 T\ = Q. (16) 

Here V and T arc the hessian matrixes of the second partial derivatives of the potential and kinetic energy respectively, 
with respect translational and rotational degrees of freedom of the three rings. 

In Tab. [l] are reported the values of the coefficients {c^} determined by the minimization together with the values for 
{fefe} derived from Eqs. (p"3|). The corresponding set of frequencies is shown in Tab. ^| where ujmd are the frequencies 
derived from a computer simulation on a system of isolated molecules at low temperature. These frequencies have 
been identified via the peaks of the spectrum of the velocity autocorrelation function. 

Other model details are reported on ref. [^(j . 



IV. COMPUTATIONAL DETAILS 



We have studied a system composed by 108 molecules (324 rings, 1944 interaction centers); the sample is large 
enough to neglect finite size effects on the investigated properties with reasonable computation times. The values 
of the parameters entering the site-site Lennard- Jones interactions were determined by preliminary simulations: the 
value of e were firstly determined by comparing the computed and the experimental |2l| ] self-diffusion coefficients 
versus temperature in the range 380 K < T < 440 K; the value of a was estimated by tuning the static structure 
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factor so to place the first maximum in the right position [E2|. Successive iterations led to e/ks = 14 K and a = 0.4 
nm. To speed up the calculation of the intermolecular potential and to assure all the torques to be estimated in 
a consistent way, the cut off for L-J interactions is applied between rings centers, i.e. the sites are considered not 
interacting when they pertain to rings whose center distance is larger than r c ~ 3(tr + L a ) = 1.6 nm. 

At this stage a remark about the rotational potential expressed in Eq. (|^) is needed: such potential presents a 
serious drawback since the intrinsic ambiguity due to the parity of the scalar products involved can force the side 
rings in the wrong positions. This problem can be bypassed introducing a new term in the intramolecular potential 
that gives zero contribution to the energy when the molecule is close to the equilibrium position. If we define the 
product w = (lx • n2){h ■ hs), and we put c-j = IOcq ( the actual value of c-j is irrelevant, it must only be able to force 
the rings in the right way ) we can write 

y _ f c 7 w 2 if w > , 17 > 
A 1 otherwise ^ ' 

This term has been activated only during the preparation of the initial configuration (i. e. at high temperature); 
during the thermalized evolution of the system we expect molecules do not drive too much away from their equilibrium 
configuration and then this term to be always equal to zero. 

To integrate the equations of motion we have treated each ring as a separate rigid body, identified by the position of 
its center of mass and by its orientation expressed in terms of quaternions |2[| . The standard Verlet leap-frog 
algorithm p3fl has been used to integrate the translational motion while, for the most difficult orientational part, the 
refined algorithm due to Ruocco and Sampoli jgj] has been employed; such choices allow a very stable integration 
with a relatively long time step. 

The rotational dynamical problem can be written as 



dJ, 



dt 
dt 



(18) 

= M{{q ii }).J ii 



where is the angular momentum of ring £ in the molecule i, f ^ the torque acting on it and M.^ is the inverse of 
the inertia tensor in quaternion coordinates. 

The expression of the torques is simplified as the rotational part V ro t of the intramolecular potential energy has 
been expressed in terms of suitable scalar products written in the general form s = vi/\vi\ ■ ^2/ 1 ^2 1- Therefore the 
value of the torques t Vi and t V2 associated to the degree of freedom corresponding to the angle arccos(s) can then be 
evaluated by ]19| 

dV rot ,_ _ , . , 

r vi = -t V2 = x "2) (19) 

The leap-frog algorithm for the rotational motion reported can be written in the form |24j 

J ie (t + t') =J ii U-^j+ fit (*) + *') (20) 
with < t' < At 

ftc(* + T) =5ie(t)+ / dt'M(q i£ (t + t'))-J it :(t + t') (21) 
Jo 

with < t < At 

J,;? ( f + t) = J ^ (* ~ t) + fl ^ t)At (22) 

The crucial point is that the dependence of the matrix M.^ on the angular variables implies the need of using a time 
step At' to perform the numerical integration appearing in Eq. (^l|), smaller then At used for the CMs integration. 
In turn, the value of At is limited by the highest vibrational frequency of about 450 cto -1 . In our case, the chosen 
values of At = 2 fs and At' = At/5 are found to be sufficiently small to reduce the fluctuations of the total energy to 
a negligible fraction of the kinetic energy. It is important to note that the use of the smaller At' does not increase 
significantly the CPU computational time since the time consuming par t in each MD step is the calculation of the 
forces and torques which are kept fixed during the integration of Eq. (EOh . 
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We considered a wide temperature range spanning the liquid phase and reaching the region close to T c as shown in 
Tab. |§| in which the whole set of simulation times is reported. During the different temperature runs, the size of the 
cubic box has been rescaled in order to keep the system at the experimental density which, for T > T g , can be fitted 
by the polynomium |2q] 

p(T) = 1.2983 - 7.00 • 10~ 4 • T - 1.23 • 10~ 7 • T 2 (23) 
where p is in g/ 'cm 3 and T in K. At each temperature T we have organized the simulations following this scheme: 

• The system was coupled for a time t resci chosen in a somehow arbitrary way, to a stochastic heat bath, i.e., the 
velocities of the rings were replaced following a logarithmic pattern with the velocities drawn from a Boltzmann 
distribution corresponding to such temperature; 

• At this point the system was at the desired temperature and we let it to perform a microcanonical time evolution 
(constant energy) for a period t te rm comparable with the experimental structural relaxation time r a at the same 
temperature; in such a way we expect every slow degree of freedom of the system to be correctly thermalized 
and we control that there was no drift in temperature and the degree of energy conservation (fluctuations are 
always less than 1% of the kinetic energy). 

• We considered the final system configuration obtained in this way as a good equilibrium starting state for a 
molecular dynamics trajectory. At each temperature we perform three different runs: the first one 20 ps long 
with a 4 x 10 -3 ps configuration saving time has been used to compute the small time dynamical behavior of the 
system; a second one 640 ps long with a 4 x 10~ 2 ps saving time has been used for some intermediate frequency 
analysis. The last one, with length t proc i and saving time t save dependent on temperature, has been used to 
calculate static quantities and the long time behavior of the system. 

All the calculations have been performed on a cluster of four a-CPU with a frequency of 500 MHz; every nanosecond 
of simulated dynamics needed approximately 24 hours of CPU-time. 



V. RECALL OF MAIN RESULTS FROM THE MODE-COUPLING THEORY 



A great improvement to our understanding of the glassy state of matter has come from the extension of the 
theoretical building of the Mode Coupling Theory (MCT) developed for the equilibrium description of the 

dynamics of simple, i.e. monatomic, liquid to the study of the glassy state. Although such theory is the only 
microscopic approach to the glass transition leading to many predictions on the experimental data, it is still at the 
center of a strong debate and some questions stay open. Infact, even if the real range of validity of MCT for the study 
of molecular liquids has been cleared in the last years (see, among others, refs. p^-p8j), some experimental results, like 
the presence of the so-called knee characterizing the low- frequency behavior of the light scattering susceptibility [£9 30 
or the presence of a cusp in the non-ergodicity parameter pjj ], seem to contradict fundamental predictions of the 
idealized version of MCT. 

In this section we sum up the main predictions of the so called ideal MCT where it is hypothesized a complete 
dynamical freezing and the so called "thermally activated hopping" processes are neglected; such predictions will be 
compared with our simulation data. In the ideal MCT the glass formation is interpreted as a dynamical transition 
from an ergodic to a non-ergodic behavior at a cross-over temperature T c . MCT provides a self-consistent dynamical 
treatment Q for the density correlation function of an isotropic system 

F(q,t) = ±(5p*(t)8p* q (0)) (24) 

where N is the number of the particles, Sp = p— < p >, 5p q {t) = Yl^—i exp(iq-r~i(t)) and r~i(i) is the position of particle 
i at time t. MCT proposes a particular ansatz for the memory kernel in the related integro-differential generalized 
Langevin equation, such kernel is coupling non-linearly the density fluctuations with one another. If the coupling 
increases upon lowering the temperature, the resulting dynamical feedback leads to a progressive slowing down of the 
density fluctuations until they become completely frozen at the critical temperature T c . The ideal MCT describes 
the behavior as much as the temperature approaches T c , i.e. the parameter a — (T c — T)/T c is small (however real 
comparisons have to be made for a not too small, in contrast to the case of scaling laws in phase transitions) . For 
temperature T > T c , F(q,t) is characterized by two step decays taking place at different time scales and the theory 
gives specific predictions for such different time regions. 
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The first one, the [3-process region, is centered around a time r CT which is predicted to scale like r CT oc \T — T c \^l 2a 
with < a < 0.5 and to be bounded in the interval tq << T a << r a where To is the time scale of the microscopic 
dynamics and t q is the structural rearrangement time scale. In this region the factorization property holds, in the 
sense that the density correlation function can be written as 



F(q,t) = f{q) + h{qU\u\G ± {t/T a ) (25) 

where f(q) is the non-ergodicity parameter (i.e. Debye-Waller factor for collective correlators or Mossbauer-Lamb 
factor for single-particle correlators), h(q) is an amplitude independent of temperature and time and the ± in G± 
corresponds to time larger or smaller with respect to t„. So, the time dependence of the correlation functions is 
all embedded in the q-independent function G±, namely spatial and temporal correlations result to be completely 
independent. G±(t) is asymptotically expressed by two power laws, respectively the critical decay and the von 
Schweidler law Q , characterized by the temperature and momentum independent exponents a and b 

G ± (-)=\ W f IjK ^«*«^ (26) 

VW \~{t/T a ) b T a «t «T a ; v ' 

here a is the same exponent of the power divergence of T a at T c and it is related to the exponent b (0 < b < 1) via 
the equation 

r 2 (i-«) = r 2 (i + b) 

r(i-2<z) r(i + 26) 1 ' 

where Y is the gamma function. 

The second time region is the so-called a-region where the second decaying step takes place. This region is connected 
to the collective structural relaxations and asymptotically the theory predicts the validity of the well known time- 
temperature superposition principle; it states that, on time scales of same order of magnitude of T a , the following 
scaling law holds at every temperature T 

F( " t)= ^(^) i (28) 

in other words, the correlation functions of any observables at different temperatures can be collapsed into a master 
curve when the time is scaled with t/r a . Moreover MCT predicts that this master curve can be fitted by a Kolrausch- 
William- Watts function {stretched exponential) 

F(q,t)~f(q) expj- (J-^ j (29) 

The a time scale r a depends on temperature trough a power law of the form 

r a cx (T ~ T c )-t (30) 
where the q-independent exponent 7 is related to the power exponents a and b of the /3-region by the relation 

11 

7=^ + ^7- 31 
la lb 

The inverse of the diffusion constant D~ 1 (T) is predicted to scale like r Q || and consequently it follows Eq. |p0|| . 

Up to now, all dynamical results reviewed arc universal in the sense that they are predicted to hold for the correlators 
of every observable with non-zero overlap with density; in particular this is true for both the one-particle and the 
collective density correlation functions. Nevertheless important differences are predicted to hold for the q-dependence 
in these two cases: in the former case f(q) and h(q) depend smoothly on q, in the latter one they oscillate respectively 
in phase and out of phase with the static structure factor S(q). Moreover, [3 a is predicted to be a smooth function of 
q in the one-particle case; at variance, it shows pronounced oscillations in phase with S(q) in the collective case. 
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VI. RESULTS 



A. THERMODYNAMICS 



In this subsection some thermodynamical time-independent results are shown like potential, kinetic, and total energy 
(see Tab. || and Fig. The interest of these results is clarified by the following argument: in computer simulations 
dealing with the glass transition it is possible to define a temperature often named T g - S i m [ p2| at which one-time 
quantities show some sort of discontinuity. Such discontinuity, whose position depends on the thermal history of the 
system, represents the thermodynamical point at which the system undergoes a glass transition on the time scale of 
the computer simulation, falling out of equilibrium. It is clear from Fig. |]that no discontinuity is present, i.e. T s _ s j TO 
of our simulations is less then the lowest temperature studied and then we have good chance to have well-thermalized 
results. 

Nevertheless, whether or not the system is in equilibrium can be checked only a-posteriori by comparing the total 
simulation time with the measured relaxation time. 

From the linearity of E(T) we deduce a specific heat c(T) constant in the temperature range investigated and equal 
to 140 ± JK~ 1 mol~ 1 ; such value must be compared with the experimental value of 341.7 JK~ 1 mol~ 1 [[53|. It is 
possible to explain the inconsistency between the two values keeping in mind that our MD value is a classical (non- 
quantic) result and, more important, we are neglecting many (~ 78) degrees of freedom concerning the deformations 
of the phenyl rings. 



B. STRUCTURE 



In general the static structure of a fluid is well described by the pair distribution function 

V_ 



« = ^EE^-I%l)>. (32) 



In computer simulations p3|| , we can identify the distances \fij\ with different quantities. 

In Fig. H we report some <?(r)'s at T = 300 K where we have considered as \fij | the distances between the Carbon atoms 
belonging to different rings (A) and between the center of masses of rings (B) and molecules (C); both total (solid 
line) and intcrmolccular (dashed line) contributions are shown in order to separate the internal molecular structure 
and the mean structural organization of the whole bulk sample. In Fig. || (B) a two peak structure is present: the 
first sharp peak is placed at r ~ 0.42 nm corresponding to the mean distance between rings belonging to the same 
molecule; the second one, of intermolecular origin, is placed at r ~ 0.6 nm. It is worth noting that such distance is 
less than the greatest intramolecular C-C distance (~ 0.7 nm ). Moreover the molecular centers of mass g(r) also 
shows a large value on distances less than 0.7 nm giving the evidence of a strong packing of the molecules. All these 
features appear to be approximately temperature independent. Such packing depends strongly on the orientational 
internal configuration of the molecules namely on the positions of the two side rings with respect to the parent one; 
the computation of the probability distribution of the scalar products among the versors l^, rh^ and n^, introduced 
in sect. O, is somehow instructive in this sense. 

In Fig. q the distribution functions for the quantities l\ ■ h 2 ^, h% -h^, fix ■ n 2j 3 are shown. The first two distributions 
are practically temperature independent and give us only informations on the correctness of the simulated geometry: 
they are sharply peaked on the the correct equilibrium positions of about 0.71 and 0.69 respectively. At variance with 
the distribution of x — h-x ■ h-$ and of hi ■ U2,3 that are symmetric around x = 0, the distribution of li ■ 77.2,3 does not 
present the symmetric peak on negative values so that we can argue that the auxiliary term Va worked correctly. The 
most interesting distribution is the third one in which the peak intensity (the peak is correctly placed at «o = 0.59) 
is higher the lower the temperature, as shown in Fig. [?], indicating therefore that the correspondent degree of freedom 
(the in-phase motion of hi ■ ^2,3) is more and more frozen on its equilibrium value with decreasing the temperature. 

We have seen that in the isolated molecule the rotational motion of the two side rings can be separated in two 
contributions: an out-of-phase harmonic libration and an in-phase pivoting around the bonds which permits rings 
to cross from one equilibrium position to the other degenerate one. It is clear from the structure of the distribution 
function in proximity of hi ■ 722,3 = shown in Fig. || that the time needed for the transition from a minimum to the 
other one will be longer lowering the temperature; moreover Fig. ^| can be considered as a restatement of the energy 
map shown in Fig. 0, being the intensity of the maximum in zero a measure of the transition probability between 
the two minima. Such fenomenology will be clarified in future communications where we will study the relaxation 
processes associated to the angular degrees of freedom. 
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The space Fourier transform of g(r) is the static structure factor. In a poliatomic system this quantity is defined as 

%)cx^£kMe^- f ^>; (33) 



N ■ 



where the coefficients bi are the scattering lengths in principle different for each species involved 



The S(q) has been determined experimentally for OTP by neutron scattering 22 39], and the following main features 
have been observed: 

i) in contrast to atomic systems its main peak is splitted in two sub-peaks placed around 14 imr 1 and 19 nm -1 

ii) in the q — > region, lowering the temperature, a reduction of scattering intensity is observed due to the decrease 
of the isothermal compressibility \t (xt oc S(q = 0)) 

Hi) increasing the density, a slight shift of the peak position to higher q values is observed 

iv) on decreasing temperature, the height of the peak around 19 nm -1 increases while the intensity of the peak at 14 
nm -1 remains nearly unaffected except for a slight reduction mostly connected to the decrease in xt- 

In Fig. ^| we show our results for the structure factors calculated assuming as scattering centers the molecules and 
the rings centers of mass with £>j = 1; every point is an average on all the independent Miller indices corresponding 
to the given q. 

It is much more interesting to make a comparison among the MD results and the experimental data obtained by 
neutron (Fig. [l(] (A)) and to test what is expected for X-Ray scattering (Fig. [llj (B)). 

In evaluating S(q) by computer simulation for a comparison with neutron data, we have to take in account the 
contribution due to both Carbon and Hydrogen atoms; H atoms are not considered in our dynamics, nevertheless it 
is possible to place them in fixed positions on the line extending from the center of the ring trough a carbon atom at 
fixed C — H distance computed to be al c-H = 0.107 nm. In this case we would have to consider different scattering 
lengths for the two species, bu and be', nevertheless they are both positive and about of the same magnitude so that 
the product bibj in Eq. ( p3[ ) is an ineffective constant. 

In Fig. [l(] (A) the calculated S(q) at T = 300 K is shown and compared with the data of ref. 0atT= 324 K; 
in this paper the authors show their results in terms of the coherent scattering cross section (dajdfl) measured in 
m 2 which is proportional to our S(q). In order to compare the two results we renormalized the experimental data in 
such a way the values of the two curves coincide at large q. The high-q region of the calculated S(q) appears to be 
in excellent agreement with the experiment but no double peak structure is present at low momenta. In particular 
the MD calculated first peak presents a small bump at about 18 nm -1 ; this is better seen in Fig. [ll] where we show 
the small-g part of S(q) calculated at T = 280 K together with the error bars estimated by means of the statistical 
fluctuation of the data. The noise cannot allow us to determine the correct structure of the main peak. It is worth 
noting, moreover, that at this low temperature the characteristic relaxation time is of order 1 ns, so that, considering 
a simulation run 10 ns long, we have only about 10 really independent system configurations. 

In order to calculate the simulated S(q) as it is expected by X-ray scattering, we consider only the Carbon atoms; 
also in this case no double-peak structure is observed in the data but a clear prepeak appears at a q-value less then 
the q of the first maximum, being the high-q behavior similar to the neutron case, as shown in Fig. Il0| (B). 



C. SELF DIFFUSION COEFFICIENT 



An important quantity to consider in the study of the dynamics of our system at a microscopic level is the mean 
squared displacement (MSD) defined as 

(r 2 (*)) = ^£<|iM*)-^(°)| 2 > (34) 

where Ri^(t) is the position of the center of mass of the ring £ in the molecule i at time t; from the MSD is possible 
to determine the self- diffusion coefficient D(T) via the Einstein relation 

£>=Hm ^<r 2 (t)) (35) 

The temperature dependence of the MSD is shown in Fig. [l^; each curve follows the usual cage-effect scenario. 
At small time (less than 0.2 ps) they present the t 2 behavior corresponding to the ballistic motion; at long time the 
diffusive linear time dependence of Eq. (|3^) is found. At intermediate times a small region is present where MSD 
stays almost constant and whose duration increases with decreasing temperature; on these time scales molecules are 
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trapped in cages builded up by their neighbors, and they can only vibrate in these limited region, the length of the 
plateau being a measure of the mean life-time of the cages. 

The calculated values of the self-diffusion coefficient are shown in Tab. |5| and plotted in Fig [l3| (open circles) as a 
function of temperature, together with the power-law temperature dependence (solid line) predicted by the MCT 

D-\T) cx (T-T c )- 7 . (36) 

A three parameters fit to these data has been performed obtaining the following values 

T^ D) = 278 ± 3 (37) 
7 (I5) = 1.8 ±0.1 (38) 

In the same figure we also show the experimental data (full squares) |2l],|3(| that are well represented by Eq. ( |3^ ) 
(dashed line) with the values 7 = 2.3 ± 0.1 and T c — 292 ± 2 K. 

It is clear from these values and from Fig. [l3] that a discrepancy is present among the lower temperatures diffusive 
behavior of the simulated and real system respectively; this is most likely due to the fact that we have tuned the 
value of the L-J potential parameters e and a in order to reproduce the high-T diffusion properties of the real system. 
However it is worth noting that it is possible to reproduce quite well the experimental results on the whole investigated 
temperature range shifting the molecular dynamics points at temperatures ~ 20 K above their true values. In other 
words we have to assume that our actual thermodynamic point is shifted with respect to the real one; from now on, 
whenever we will compare our molecular dynamics results with the experimental ones, our calculated points will be 
shifted of 20 K above the measured temperature and the competing temperatures will be indicated as T c . On this 
grounds from the previous study of the self diffusion properties of our model we obtain Tc D%> = 298 ± 3 to be compared 
with the experimental value T c = 290 ± 5. 

A different way to determine the parameters entering in the power law reported in Eq. [36] is possible, even if not 
independent from the previous one; it is based on the study of the non-gaussian parameter 0:2 (t) defined as pftpq]; 

9 (r 4 (t)) 

where the mean square displacement (r 2 (t)) and (r 4 (t)) are, respectively, the second and forth momenta of the Van 
Hove self- correlation function 

G ^ = i£(*(r--R«(t) + £<e(0))>. (40) 

»£ 

The parameter c*2(i) quantifies the degree of non-Gaussianity of G s (f,t) in space as a function of time and it is 
normalized in such a way that, if G s (f, t ) was a gaussian function in space at a given time t, we would have 02 (t) = 0. 
The time dependence of «2 (t) at all temperatures investigated is shown in Fig. |lj. We are not interested here in the 
specific time dependence of such function but only in the fact that t a „ , the position of the maximum of a.i (t), has 
the power dependence on T similar to that of Eq. ( |36j ) |37j (see Fig. [l||) . A fit to this quantity performed in the same 
way as before gives us the values 

T c tQ = 300±14 (41) 
7*° =1.4 ±0.3 (42) 

compatible with the values determined by the temperature dependence of D, even if the error bars are larger in this 



D. SINGLE PARTICLE DYNAMICS 



Comparisons of the coherent (collective) and incoherent (self) density fluctuations dynamics data measured by differ- 
ent techniques (neutron time-of-flight and backscattering spectroscopy, photon correlation spectroscopy, depolarized 
Raman and Rayleiah-Brillouin light scattering) with the main predictions of MCT have been reported in literature 
with great details ||,[3!^-|^| . In this section we will study the single particle density fluctuation dynamics of our model 
and we will compare our results with the experimental results mainly contained in refs. pj,|39[] and with the MCT 
predictions. 
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The single particle dynamics of the model is embedded in the incoherent self-intermediate scattering function defined 

as 

Fs ^ t ) = l^ e -*[%W-«*e(0)]) (43) 

where, again, (t) is the position of the center of mass of the ring £ in the molecule i at time t. At every temperature 
considered two sets of configurations, produced with two different storing time as described in section [V, have been 
used to reconstruct the whole curve. We considered the T-dependence of F s (q,t) at the two momentum values 
q = 14, 19 nm^ 1 corresponding to the first and second peaks of the static structure factor, averaging on values of q 
falling in the interval q ± Aq with Aq — 0.2 ran -1 . Finally, we spanned at T = 300 K the whole interesting g-space 
in the interval q — 2 -j- 30 nm _1 (averaging on the values of q falling in the same interval 2Aq wide). 



In Fig. 16 we show F s (q,t) for nearly all temperatures investigated at q = q max — 14 nm _1 ; all the curves decay 
to zero i.e. the length of all the simulations allows the fluctuations to become completely uncorrected. We are in a 
"good" thcrmodynamical equilibrium at every temperature, at least on the space scales corresponding to the inverse 

of Qmax • 

At temperatures lower then T = 330 K the relaxation follows clearly the predicted two step pattern: on microscopic 
time-scales the correlation is quadratic in time, this time scale being the one on which the intramolecular vibrations 
happen; on intermediate time scales we observe the formation of a plateau, whose height is the non-ergodicity param- 
eter f(q) and whose length in time is comparable to the one of the plateau in the MSD (r 2 (t)). On long time scales 
we observe the structural relaxation in the form of a stretched exponential. At the highest temperatures no double 
pattern is visible anymore and only a nearly exponential relaxation can be recognized. A stretched exponential fit (see 
Eq. (|2"g|)) on the structural time scale (a-process) gives us the temperature dependence of the three free parameters 
[3 a ,T a and f(q). 

The parameter (3 a (circles) is shown in Fig. |lj; it appears to be nearly T-independent for temperature lower than 
T = 400 K and its mean value [3 a ~ 0.8 (dashed line) has to be compared with the experimental value [3 a = 0.6. 
For temperature in the higher region it tends toward the value (3 a = 1 (errors are clearly much more greater); such 
behavior is due to the fact that in this temperature region it is no longer possible to sharply separate the long-time 
relaxation region from the microscopic short-time one. The study of the temperature dependence of the non-ergodicity 
parameter f(q) in the interesting region is not possible due to our limited temperature range which do not permit 
to observe the expected low-temperature (T < T g ) harmonic Debye-like behavior, and the onset of the anomalous 
decrease of f(q) with increasing T for T g < T < T c . Our data suggest us T c < 283 K and the mean value f(q) ~ 0.7 
(dot-dashed line) agrees with the experimental value determined at T = 290 K shown in Fig. ^l|. 

It is worth to test the power law temperature dependence Eq. ( pO] ) for the relaxation time r Q ; the calculated 
relaxation times (circles) shifted of 20 K with respect to the measured temperatures, as explained above, are plotted 
in Fig. |l8] together with the experimental (full triangles) shear viscosity Tj s data (r] s is expected to be proportional to 
T a ) of ref. p5| and the theoretical fitted curve (solid line) of parameters 

T> = 296 ± 7 (44) 
7 s = 2.0 ±0.4 (45) 

to be compared with the experimental results of ref. (|] T c = 290 ± 5 K and 7 = 2.55. These values are compatible, 
within the statistical error, with the values calculated from the diffusion data so that we can conclude that the diffusive 
behavior and the self-dynamics of our model follow the same critical power law with T c ~ 297 K and 7 ~ 1.8. 

Also the values of T a for q = 19 nm _1 (squares) corresponding to the second peak of the static structure function 
are reported with the theoretical curve. A fit has been performed (dashed line) only on the prefactor keeping fixed 
the values of the other two parameters in order to show that these data also are compatible with the same power law. 
The crucial observation here is that the values of the two parameters T c and 7 are effectively q-independent and they 
can be considered universal for our model, as predicted by the MCT. 

The relaxation time r Q can be also used to test the time-temperature superposition principle Eq. (psj). In Fig. ^ 
the curves are shown in function of the rescaled time t/r a and it is clearly seen that all the curves tend to collapse 
on the same master curve as predicted by the theory. 

We now quantify the q-dependence of the self dynamics long-time behavior of the system at T = 300 K. In Fig. |2^ 
are reported the curves F s (q,t) for values of q = 2n nm _1 with n = 3, . .., 15; the choice of the temperature value 
T = 300 K has been due to the need of "well-thermalized" results in a large range of q. Also in these data is well 
defined the two step behavior and we can calculate the long time stretched exponential fit parameters; the resulting 



values are shown in Fig. |2l| and Fig. |22j In Fig. |21]the values of f3 a (left side) and f(q) (right side) are shown. 

P a (open circles) appears to be a smooth function of q and it tends, for large values of q, to the experimental (full 
circles) evaluated value of 0.6. Such behavior is quite general (see, for instance, ref. |38|) and can be easily explained 
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by the following argument |38|: for large values of q, corresponding to length scales of the same order of magnitude 
of the cages dimension, the dynamics becomes slower and slower approaching the cage dynamics described trough 
the von Schweidler exponent b. At variance, in the opposite limit of small q, we consider a diffusive dynamics on 
large distances; at such length scales the decay of the self-density fluctuations is of the usual purely exponential form 
exp (—Dq 2 t) corresponding to j3 a = 1 (see Fig. ^). At this stage, however, we have not a reasonable explanation of 
the disagreement with the experimental data. 

In the right side of Fig. [2l] the q-dependence of the non-ergodicity parameter f(q) is also shown as calculated from 
the short time limit of the a-process (squares), from the long-time part (triangles) of the /3-region (as we will see 
below) and from the experimental data (full squares) ||] ; it seems clear a good agreement between our values and the 
experimental results. The data appear to be monotonic decreasing as increasing q and this dependence is expected 
to be approximately gaussian; in Fig. ^l] a gaussian fit (solid line) in the form exp(— q 2 /2a 2 ) with a ~ 19 nm _1 to 
molecular dynamics a-region data is also shown. It is clear that the g-range considered here is too limited to really 
decide on the validity of this functional form (a linear approximation would work well too), a good estimate of the 
error bars lacking in this case. 

In Fig. ^2] the g-dependence of r" 1 is shown (circles) together with the experimental data §|(full circles). Molecular 
dynamics points have been rescaled by a factor T^f D (T = 280K)/t£ 1d (T — 300K) ~ 6.5 to take into account the fact, 
as discussed above, that our system temperature is 20 K higher the real one. The correct square law behavior at low-g 
r Q T 1 (g) ~ 6Dq 2 (see Eq. [l7]) is also shown as a solid line; here D is the self diffusion coefficient and 6D — 20 .4 x 10~ 5 
nm 2 /ps. 

Finally, in Fig. E3, two products of some calculated quantity, expected to be constant, are shown; on the left side 
the statement D~^T) cx r a (T) is proven for q = q max (good at highest temperature). On the right side we show 
the product r a q 2 D evaluated at T = 300 K which is nearly constant up to q ~ 18 as expected, being this value 
approximately the crossover point among the correct quadratic behavior and the asymptotically linear regime as we 
found above. 

Let us now probe the MCT conclusions about the /3-region which is predicted to follow the power laws of Eq. (p6|): 
we show here the results at fixed temperature T — 300 K for values of q — 6 + n nm -1 with n — 0, . . . , 17. 
We fit all the curves by two power functions on the time ranges 0\ and 62'- 

F ,„ .x / f q + c 1 t- a ted! 
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where 6\ = [0.15 : 2] (ps) and O2 — [3 : 20] (ps); some selected fits are shown in Fig. |2J and they seems to work quite 
well. 

Some observations arc needed on the following analysis: a great uncertainty stems from the choice of the fit range 
(i.e. from the choice of #1 and 62) due to the consideration of a cross-over region between two process not sharply 
separated in time; moreover it cannot be excluded an analogous problem between the microscopic region and the 
critical decay region characterized by the exponent a. Such difficulty implies a great uncertainty on the determination 
of f(q) which is supposed to be considered as the long-time limit of the (3- process and the short-time behavior of the 
a-relaxation. At this point it is clear that, lacking a careful error analysis on the data points, such fits can only state 
that a parameterization of the data in the form of Eq. (^6|) is possible which is consistent with the theory predictions 
f39f , the current values of the determined parameters being considered only from a qualitative point of view. 

Nevertheless the values of such fitting parameters, shown in Fig. |2^ and Fig. p^, appear to be in good agreement 
with the values obtained from the experimental data. In Fig. ^5] are shown (left)the values of the power exponents 
a (squares) and b (circles) and of the exponent 7 (triangles) calculated by means of Eq. (31) (right side); the mean 
values are 0.3, 0.5, 2.6 respectively, to be compared with the experimental determined values a ~ 0.31 (dot-dashed 
line), b ~ 0.52 (dashed line), 7 ~ 2.55 (solid line) |]. 

It is clear from these results that one of the main prediction of MCT, namely the g-independence of a and b 
is verified in the limit of the error fluctuations. Moreover it is important to note that the parameter 7, given by 
Eq. (|3l|), remains constant, as expected, being its mean value 2.6 clearly compatible with the value 2.55 determined 
from experimental data ||; nevertheless this value overestimates the value 7 = 2.0 ± 0.4 at q = 14 nm _1 previously 
determined by the fit to the a-region. Furthermore b(q) is always less then the determined value of the large-g value 
of the stretching parameter (3 a = 0.6 (see Fig. |2l]) verifying an other MCT prediction, namely < b < (3 a . 

From Eqs. ( |25| ) ( |2^ ) ( pH ) the two parameters c\(q) and €2(5) result to be proportional to h(q), the proportionality 
constants being dependent on a, r CT , a and b. From T c ~ 280 K we have \fa ~ 0.3 while we choose as a good estimate 
of T a the intersection point of the two power laws of Fig. |24| obtaining r a ~ 2 ps; if we put a — 0.31 and b = 0.52 we 
finally obtain the factors 2.7 and 2.3 for C\ and C2 respectively. 

Unfortunately these value are not able to correctly rescale our data on the experimental results, the correct values 
being 0.7 and 4 as shown in Fig. E6c this result could be expected considering the great uncertainties on the parameter 
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values used for the estimate. Nevertheless, simulated data agree quite well with the experimental points at low-q 
presenting a strong bending toward a constant value in the region q > 16 nm . 

To complete the picture of the self-motion in our model, we test the validity of the gaussian approximation to 
F s (q, t) in the limit of small momentum q. The first order term of the expansion of F 8 (q, t) in powers of q 2 gives |34| l 

F s M)cexpj-^(r 2 (;))j. (47) 

In Fig. |3 some curves at T = 330 K and different values of q are shown together with the corresponding approxima- 
tions; such approximation seems to work quite good and it becomes worse on increasing q as expected. 

At the end, we show in Fig. ^8] all the time scales related to centers of mass motion investigated up to now as a 
function of temperature. Full circles and squares indicate respectively the experimental structural relaxation time 
tYxp(T), following the Vogel-Fulcher law, and the experimental inverse self-diffusion coefficient D^} xp {T) multiplied 
by a factor 5 x 10~ 5 cm 2 in order to superimpose to T E xp{T). The open symbols are used to represent the molecular 
dynamics results: t^ d (diamonds) is the relaxation time of the one-particle dynamics at q = 14 nm" 1 multiplied by 
a scale factor 1.5 and D^ D {T) (triangles up) is the inverse of the diffusion coefficient rescaled by the same factor 
5 x 10~ 5 cm 2 we used for the experimental points. All the molecular dynamics points have been shifted of 20 K with 
respect to the measured temperatures; it is quite clear that both experimental and molecular dynamics data points 
collapse on a well defined master curve. Our model is, at least, a good model for centers of mass dynamics of OTP. 



VII. CONCLUSIONS 



In this paper we have introduced a new interaction potential model capable to describe the intramolecular dynamics 
of the fragile glass former OTP; such model appears to be much more efficient with respect to the ones introduced so 
far in the sense that it represents a much more better compromise between the resulting computing needing and its 
capability to mimic all the complexity of the dynamical behavior of the real system. 

It takes into account not only the translational and rotational dynamics of the molecules as a whole but also the 
stretching along the molecular bonds, the tilt of the bonds, the rotations of the side rings with respect to the parent 
ring. It is tuned in such a way to reproduce the isolated molecule vibrational spectrum. In this way, most likely, we 
have introduced the degrees of freedom whose interplay causes the complex dynamical behavior of the real system. 

We have, then, presented the results of molecular dynamics computer simulations of such model; we mainly studied 
the static structure of a bulk sample, the self diffusion properties and the self part of the density-density correlation 
functions. 

The static structure factor simulated in such a way is compared with the experimental measures and shows a good 
agreement with the neutron scattering data except for the very low momenta behavior (due, probably, to the finite size 
of our system). Moreover we have no evidence of the splitting of the main peak in two subpeak placed at q = 14, 19 
nm , only the first one being clearly visible. This luck may be due mainly to the temperature range investigated 
(the intensity of the second sub-peak increases with lowering temperature). 

The self-diffusion properties of the system have been investigated trough the mean squared displacement and the 
self diffusion coefficient temperature behavior; comparisons with experimental self-diffusion data give a very good 
agreement, showing the evidence of compatible critical dynamics behavior in approaching the instability temperature 
T c of MCT which is here find to be T c — 278 ± 3 to be compared with the experimental value T c = 290 ± 5 determined 
by a MCT analysis of the dynamics of the density fluctuations, a discrepancy that is most to likely ascribed to 
the intermolecular LJ potential parameters (a, e) that have been tuned in the temperature region close to T = 300 
K. Moreover we considered the critical temperature dependence of the so-called non-gaussianity parameter a^if) 
obtaining compatible values for the power law parameters. 

The self dynamics of the density fluctuations has been studied in great detail on the whole accessible time window, 
spanning the range from time scale of the order of few fs to times of order of some ns; moreover its dependence on 
temperature and momenta has been investigated. All the correlation curves calculated show the typical two step 
behavior predicted by MCT, the first one on short time being associated to the so-called "microscopic" processes, i.e. 
the vibrational motion of molecules in the cage built up by their neighbors; the second one being associated to the 
a-process which controls the structural rearrangements of molecules on long time scale. 

The critical dynamics on a time scale approaching the correspondent T c is in good agreement with experimen- 
tal findings; indeed the estimate values for our model of the exponent 7 = 2.0 ± 0.4 must be compared with the 
experimental value 7 — 2.5. 
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The q-dependence at T = 300 K in the momentum range q = 6 — 30 nm _1 has been analyzed in terms of a stretched 
exponential fit; the values of the determined parameters are in good agreement with the ones calculated by fitting the 
experimental points and with the MCT expected behavior. 

To summarize, in the present work we have shown that our model for OTP fluid is mimicking rather well the center 
of mass dynamical features of the real system, giving results in most cases fully compatible with the experimental 
findings. It is clear that its skill to help us in the understanding of the most exotic dynamic features of the real 
systems and, in particular, of the relevance of the internal degrees of freedom on the translational dynamics, has 
not been cleared in this paper; problems like the collective density fluctuation behavior, the rotational dynamics, the 
origin of the unusual fast relaxational dynamics and many others, will be addressed in future works p0| . 

ACKNOWLEDGEMENTS 

The authors wish to thank W. Gotze for a critical reading of the manuscript, L. Bernardini and G. Giuliani of the 
"Parco Tecnologico d'Abruzzo" for technical support, G. Monaco and F. Sciortino for some useful discussions. 



[1] C. A. Angell, Science 267, 1924 (1995). 

[2] W. Gotze, in Liquids, Freezing and the Glass Transition, edited by J. P. Hansen, D. Levesque and J. Zinn- Justin (North- 
Holland, Amsterdam, 1991); W. Gotze and L. Sjorgen, Rep. Prog. Phys. 55, 241 (1992); W. Gotze, J. Phys.: Condens. 
Matter, 11, Al (1999). 

[3] R. Schilling, in Disorder Effects on Relaxational Processes edited by A. Richert and A. Blumen (Springer Verlag, 1994); 

W. Kob, in Experimental and Theoretical Approaches to Supercool ed Liquids: Advances and Novel Applications edited by 

J. Fourkas et al. (ACS Books, Washington, 1997). 
[4] K. 1. Ngai and K. Y. Tsang, Phys. Rev. E 60, 4511 (1999). 
[5] M. Mezard and G. Parisi, Phys. Rev. L ett. 82, 747 (1999) . 



[6] F. Sciortino, W. Kob and P. Tartaglia, :ond-mat 9911062, (1999) 



[7] G. Monaco, D. Fioretto, C. Masciovecchio, G. Ruocco and F. Sette, Phys. Rev. Lett. 82, 1776 (1999). 

[8] W. Petry, E. Bartsch, F. Fujara, M. Kiebel, H. Sillescu and B. Farago, Z. Phys. B 83, 175 (1991). 

[9] G. M. Brown and H. A. Levy, Acta Crystallogr. B 35, 785 (1979). 
[10] S. Aikawa, Y. Maruyama, Y. Ohashi and Y. Sasada, Acta Crystallogr. B 34, 2901 (1978). 
[11] A. Criado, F. J. Bermejo, A. de Andres and J. L. Martinez, Mol. Phys. 82, 787 (1994). 
[12] J. J. Ou and S. H. Chen, J. Comp. Chem. 19, 86 (1998). 

[13] G. Wahnstrom and L. J. Lewis, Physica A 201, 150 (1993); L. J. Lewis and G. Wahnstrom, Solid State Comm. 86, 295 
(1993); L. J. Lewis and G. Wahnstrom, J. of Non-Crystalline Solids 172-174, 69 (1994); L. J. Lewis and G. Wahnstrom, 
Phys. Rev. E 50, 3865 (1994). 

[14] C. M. Roland, K. L. Ngai and L. J. Lewis, J. Chem. Phys. 103, 4632 (1995); C. M. Roland and K. L. Ngai, J. Chem. 

Phys. 106, 1187 (1997); F. Sciortino and P. Tartaglia, J. Phys.: Condens. Matter, 11, A261 (1999). 
[15] S. R. Kudchadkar and J. M. Wiest, J. Chem. Phys. 103, 8566 (1995). 
[16] W. R. Busing, J. Am. Chem. Soc. 104, 4829 (1982). 
[17] D. H. Wertz and N. L. Allinger, Tetrahedron 35, 3 (1979). 

[18] G. E. Forsythe, M. A. Malcolm and C. B. Moler, Computer Methods for Mathematical Computations (Prentice Hall 1976). 
[19] H. Goldstein, Classical Mechanics (Addison- Wesley Publishing Company, Inc., Reading, Massachusetts, 1965). 
[20] S. Mossa, PhD Thesis, unpublished (1999). 

[21] D. W. McCall, D. C. Douglass and D. R. Falcone, J. Chem. Phys. 50, 3839 (1969). 

[22] E. Bartsch, H. Bertagnolli, P. Chieux, A. David and H. Sillescu, Chem. Phys. 169, 373 (1993). 

[23] M. P. Allen and D. J. Tildesley, Computer Simulation of Liquids (Clarendon Press, Oxford, 1989). 

[24] G. Ruocco and M. Sampoli, Mol. Phys. 82, 875 (1994). 

[25] G. Monaco, PhD Thesys, unpublished. 

[26] L. Fabbian, A. Latz, R. Schilling, F. Sciortino, P. Tartaglia and C. Theis, Phys. Rev. E 60, 5768 (1999). 
[27] S. Kftmmerer, W. Kob and R. Schilling, Phys. Rev. E 56, 5450 (1997). 

[28] T. Franosch, M. Fuchs, W. Gotze, M. R. Mayr and A. P. Singh, Phys. Rev. E 56, 5659 (1997). 

[29] J. Wiedersich, T. Blochowicz, S. Benkhof, A. Kudlik, N. V. Surovtsev, C. Tschirwitz, V. N. Novikov and E. Rossler, J. 

Phys.: Condens. Matter, 11, A147 (1999). 
[30] J. Gapinski, W. Steffen, A. Patkowski, A. P. Sokolov, A. Kisliuk, U. Buchenau, M. Russina, F. Mezei and H. Schober, J. 

Chem. Phys. 110, 2312 (1999). 
[31] F. Mezei and M. Russina, J. Phys.: Condens. Matter, 11, A341 (1999). 



15 



[32] W. Kob and H. C. Andersen, Phys. Rev. E 51, 4626 (1995). 
[33] S. S. Chang and A. B. Bestul, J. Chem. Phys. 56, 503 (1972). 

[34] J. P. Hansen and J. R. Mc Donald, Theory of Simple Liquids (Academic Press Limited, London, 1986). 
[35] A. Tolle, H. Schober, J. Wuttke and F. Fujara, Phys. Rev. E 56, 809 (1997). 
[36] F. Fujara, B. Geil, H. Sillescu and G. Fleischer, Z. Phys. B 88, 195 (1992). 
[37] T. Odagaki and Y. Hiwatari, Phys. Rev. A 43, 1103 (1991). 

[38] F. Sciortino, P. Gallo, P. Tartaglia and S. H. Chen, Phys. Rev. E 54, 6331 (1996). 

[39] M. Kiebel, E. Bartsch, O. Debus, F. Fujara, W. Petry and H. Sillescu, Phys. Rev. B 45, 10301 (1992). 

[40] E. Bartsch, F. Fujara, J. F. Legrand, W. Petry, H. Sillescu and J. Wuttke, Phys. Rev. E 52, 738 (1995). 

[41] A. Tolle, J. Wuttke, H. Schober, O. G. Randl and F. Fujara, Eur. Phys. J. B 5, 231 (1998). 

[42] Y. Hwang and G. Q. Shen, J. Phys. 11, 1453 (1999). 

[43] W. Steffen, A. Patkowski, H. Closer, Phys. Rev. E 49, 2992 (1994). 



1G 



TABLE CAPTIONS 



1 . Values of the internal potential coefficients determined by a least square minimization of the error function, Eq. ( p.5| ) . The coefficients 
Cfc are the force constants associated to the different intramolecular potential terms; the b^, derived from the {c^} via the Eq. ([nj), 
describe the form of the in-phase rotational potential according to Eq. (|lo|). 

2. Free molecule vibrational frequencies: in the second column are reported the values luhf determined by a Harthree-Fock calculation 
of the ground state of the isolated molecule (we consider only the 12 lowest eigenvalues); in the third column are shown the frequencies 
t^MD corresponding to the different peaks of the spectrum determined from the atomic velocity autocorrelation function calculated 
by means of a preliminary molecular dynamics simulation of the isolated molecule (e = 0) at T = 1 K. 

3. Simulation runs details: at each temperature T, the system is coupled for a time t resc to a stochastic heat-bath, then it evolvs freely, 
at constant energy, for a time t tarm - At the end of this process, we consider the final configuration to be in a good equilibrium 
state and we start a trajectory t prol i ps long, saving a system configuration every t sa ve ps. 

4. Some thermodynamical results: temperatures T effectively measured, total energy Etot, total potential energy Vtot and kinetic 
energy T. 

5. Temperature dependence of the molecular dynamics self-diffusion coefficient D. 
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FIGURE CAPTIONS 



1. Molecular structure of OTP (C18-H14); it is constituted by three phenyl rings, the two side rings being attached to the parent (i.e. 
central) ring by covalent bonds. 

2. Model geometry: each phenyl ring is represented by a rigid hexagon of side L a = 0.139 nm and the equilibrium bond lenght is 
Lj, = 0.150 nm. Ci, C2, C3 represent the origins of the reference frames fixed with the rings; U2, «3, "1*1(2), %(3) are the vectors 
parallel to the ring bonds; l\ and rh\ are two versors identifying the parent ring plane; P2, P3, ^1(2), -^1(3) are positions of 
the Carbon atoms bonding together the rings; S2, S3, Si(2)> ^1(3) are four interaction sites introduced to force the rings towards 
the co-planar equilibrium condition; the symbols P\ and Pg have been introduced to identify the angles P2 — P\(2) ~ Pa an d 
F3-Pi( 3 )--Pb- 

3. Rotational energy: each point has been determined by an ab-initio Hartree-Fock calculation of the potential energy surface as a 
function of the rotational angles <f>\ and 02 (expressed in degrees) with all the other degrees of freedom fixed to their equilibrium 
values; the energy is expressed in Hartree units (1 Hartree=27.2eV). We quote by A the isolines ranging from 0.0505 to 0.0545 
mHartree with an increment Ai = 0.005 mHartree, by B the ones from 0.055 to 0.095 mHartree with A2 = 0.05 mHartree and by 
C those from 0.5 to 1 mHartree with A3 = 0.5 mHartree. This figure has to be considered only from a semi-quantitative point of 
view as explained in the text. 

4. Temperature dependence of the energies tabulated in Tab. ^, Etot (circles) Vtot (squares) and T (triangles up) together with the 
internal \ 'intra (triangles down) and the intermolecular Lennard- Jones Vi„t er (diamonds) contributions to Vtot- 

5. Pair static distribution functions at T = 300 K calculated on atoms (A), ring centers of mass (B) and molecular centers of mass (C); 
full lines represent the total contribution of both intramolecular and intermolecular distances, dashed lines only the intermolecular 
contribution. 

6. Static distribution functions of the scalar products fi2 -fi3 (triangles), li -712,3 (circles), and ni -712,3 (squares) evaluated at T = 280 
K. 

7. Intensity of the peaks of the static distribution function of fi\ • n.2,3 at T = 440, 390, 330, 280 K (from bottom to top). A blow-up 
of this figure around x = is shown in Fig. ^. 

8. Temperature dependence of the static distribution function of ni ■ 712.3 near the saddle point position hi ■ 712,3 = at T = 
440, 420, 390, 370, 350, 330, 320, 300, 280 K (from top to bottom). 

9. Static structure factors at T = 300 K calculated on on the molecular (solid line) and ring centers of mass (dashed line); each q 
point is the average over all the independent Miller indeces corresponding to it. 

10. Top: Comparison among molecular dynamics structure factor (solid line) calculated taking in account both Carbon and Hydrogen 
atoms as scattering centers and experimental structure factor (dashed line) measured by neutron scattering (from p^]). Bottom: 
Molecular dynamics structure factor (solid line) calculated taking in account only Carbon atoms; this should be the correct result 
to be compared with experimental structure factor measured by X-ray scattering. 

11. Enlargement of the low-momenta region of S(q) calculated at T = 280 K: the error bars are estimated from the fluctuation of the 
single configuration S(qYs. 

12. Temperature dependence of the mean square displacement {r 2 (t)) calculated on ring centers of mass at all temperature investigated 
except T = 440,420 K (higher temperature on top). Inset: linear scale plot of the mean square displacement at some selected 
temperatures (open symbols) together with the long time linear behavior (dashed lines). 

13. Temperature dependence of molecular dynamics (open circles) and experimental (full triangles) diffusion coefficients together with 
the power law fits in the form of Eq. |36] (solid and dashed lines respectively); we show also the MD data shifted of 20 K (open 
squares) as explained in the text. 

14. Time dependence of the non-gaussian parameter 012 it) for all temperature investigated (lower temperatures on top). 

15. Power law fit of the temperature dependence of the position t ma x of the maximum of the non-gaussian parameter. 

16. Temperature dependence of F s (q,t) calculated at q = 14 nm" 1 for all temperatures investigated except T = 410,430 K (lower 
temperatures on top). 

17. Temperature dependence of the stretching parameter f3 a (circles) and of the non-ergodicity parameter f q (squares); the horizontal 
lines indicate the mean values of (3 a (dashed line) and f q (dot-dashed line). 



18 



18. Temperature dependence of T a at q = 14, 19 nm _1 (circles and squares respectively) together with the power law fits with Tj" = 296 
K and 7 T = = 2.0 (solid and dashed lines respectively); also the experimental shear viscosity r) s oc r a data (full triangles) are reported 
(sec pj| and reference therein) multiplied by a factor 1.5 ps/Poise. Molecular dynamics results have been shifted of 20 K with 
respect to the measured temperatures, as explained in text. 

19. F s (q,t) at q = 14 nm rescaled to t/r a ] all the curves verify the time-temperature superposition principle. 

20. Q-dependence of F a (q,t) at T = 300 K for q = 2n nm" 1 with n = 3, ...,15 (from top to bottom). 

21. Q-dependence of the stretching and non-ergodicity parameters: Left: molecular dynamics (circles) and experimental (filled circles, 
from [g) values of the coefficient fj a as determined by the KWW fits; Right: experimental values (filled squares) of the non-ergodicity 
parameter Q together with the molecular dynamics results as determined by a MCT analysis of both a (squares) and /3 (triangles) 
regions and the gaussian fit (solid line) to a-region results with <r 2 = 365 nm -2 . 

22. Q-dependence of the inverse relaxation time t& ; molecular dynamics data (open circles) have been multiplied by a factor 6.5 in 
order to overlap the experimental data (full circles), as explained in the text. The solid line is the correct small-g behavior 
t~ 1 {q) ~ 6Dq 2 where 6D = 20.4 X 10~ 5 nm 2 /ps. 

23. Left: Temperature dependence of the product r a D at q = 1.4 nm - 1 ; Right: q-dependence at T = 300 K of r Q Dq 2 . These quantity 
are expected to be constant. 

24. Power laws in the /3-region for q = 8,10,12,14 nm" 1 ; critical decay of exponent a (dot-dashed lines) and von Scweidler law of 
exponent b (dashed lines) have been reported with simulation points. 

25. Q-dependence of the power laws exponents: Left: the experimental values (filled circles) for b(q) are plotted together with the 
molecular dynamics values for b(q) (open circles) and a(q) (open squares); the mean values for our results are shown by the 
dashed and dott-dashed line respectively. Right: our values for the q-dependence of the exponent 7 determined by the relation 
7 = l/2a + l/26. 

26. Q-dependence of the experimental (full circles) coefficient h of Eq. [25] and of our fitting parameters ci (triangles) C2 (squares); note 
that our results have been rescaled by a factor 0.7 and 5 respectively in order to superimpose to experimental data. 

27. Gaussian approximation to F a (q, t) at T = 330 K for q = 2, 3, 4, 6, 8, 1, 1.2 nm - 1 (from top to bottom). 

28. Master plot of the temperature dependence of all the centers of mass time scales discussed: we have used full symbols for experimental 
results and open symbols for molecular dynamics data. Molecular dynamics points collapse exactly on the master curve identified 
by the experimental data if they are multiplied by a scale factor (taking into account the momentum dependencies of the relaxation 
time t~md ant ^ correct dimensionality of the diffusion coefficients) and the corresponding temperatures are shifted of 20 K above 
the measured ones, as discussed in the text. In particular t^ [d (open diamonds) is multiplied by a factor 1.5, the self diffusion 
coefficients D~^} xp (full squares) and D~^ D (open triangles) are rescaled by a factor 5 X 10 — 5 cm 2 . 
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TABLES 
Tab. 1 



Cl 


1.36 


io- 16 


J/nm 2 


C2 


6.06 


lO" 1 ' 


J/nm 2 


C3 


2.44 


10 -iy 


J 


C 4 


2.65 


io- iV 


J/nrn 2 


C5 


4.23 


10 -iy 


J 


C6 


7.01 


io- lu 


J 


6l 


3.62 


10-* 


J 


62 


-4.11 


10 -iy 


J 


63 


6.92 


10" iy 


J 



Tab. 2 





u>HF(cm 1 ) 


uiMDicm' 1 ) 


1 


49. 


55. 


2 


54. 


85. 


3 


63. 


99. 


4 


97. 


148. 


5 


114. 


165. 


6 


140. 


172. 


7 


256. 


273. 


8 


271. 


302. 


9 


291. 


306. 


10 


358. 


376. 


11 


386. 


409. 


12 


428. 


436. 



Tab. 3 



T(K) 


tresc (ps) 




tprod{ps) 


tsave (ps) 


443 


100 


1000 


1000 


0.1 


433 


100 


1000 


1000 


0.1 


420 


100 


1000 


1000 


0.1 


410 


100 


1000 


1000 


0.1 


389 


200 


1400 


2000 


0.1 


372 


200 


1400 


2000 


0.1 


351 


200 


1400 


2000 


0.1 


331 


500 


2000 


4000 


0.2 


321 


500 


3000 


6000 


0.3 


313 


500 


3000 


6000 


0.3 


300 


500 


3000 


10000 


0.5 


294 


500 


5000 


10000 


0.5 


283 


500 


5000 


10000 


0.5 



20 



Tab. 4 



T(K) 


E to t(kJ /mol) 


T T (111 1\ 

Vtot{kJ/mol) 


T(kJ/mol) 


443 


30.00 


oil - i r\ to 

—3.15 ± 0.72 


33.12 ± 1.29 


433 


28.59 


o fro i r~\ t' i\ 

— 3.78 ± 0.69 


32.37 ± 1.26 


420 


26.88 


—4.59 ± 0.69 


31.47 ± 1.26 


410 


25.47 


—5.25 ± 0.69 


30.72 ± 1.20 


389 


22.44 


—6.69 ± 0.63 


29.13 ± 1.14 


372 


19.98 


—7.83 ± 0.60 


f^^T O ~l 1 T Aft 

27.81 ± 1.08 


351 


17.10 


-9.15 ± 0.57 


26.28 ± 1.02 


331 


14.34 


-10.47 ±0.54 


24.81 ±0.96 


321 


12.75 


-11.22 ±0.54 


24.00 ± 0.96 


313 


11.73 


-11.67 ±0.51 


23.43 ± 0.93 


300 


9.96 


-12.51 ±0.51 


22.47 ±0.90 


294 


9.15 


-12.90 ±0.48 


21.99 ± 0.90 


283 


7.59 


-13.59 ±0.48 


21.18 ±0.84 


Tab. 5 


T(K) 


10 7 x D(cm 2 /s) 


443 


106.2 


433 


89.3 


420 


81.3 


410 


66.6 


389 


49.3 


372 


37.3 


351 


27.3 


331 


12.5 


321 


9.4 


313 


6.3 


300 


3.4 


294 


0.8 


283 


0.6 
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Figure 4 
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